function [ conn, times, dat_seed, dat, time_indices, conn2 ] = osl_conn_hilbert_envelope( S )

try tmp=S.dat; catch, error('S.dat not set'); end;
try tmp=S.tfdat_seed; catch,  end; % e.g. Hilbert transformed seed data if available (do avoid it having to be recomputed)
try tmp=S.times_in; catch, error('S.times_in'); end;
try tmp=S.time_range; catch, error('S.time_range'); end;
try tmp=S.bc; catch, S.bc=0; end; % baseline correct within trial
try tmp=S.time_moving_av_win_size; catch, S.time_moving_av_win_size=0.1; end; % window to average over for downsampling (in secs)

time_indices=S.times_in>=S.time_range(1) & S.times_in<=S.time_range(2);

times(1) = mean(S.times_in(time_indices));
conn=zeros(size(S.dat,1),1);
conn2=zeros(size(S.dat,1),1);

for tri=1:size(S.dat,1), % indexes trials
    
    if(S.connectivity_seed_regress_zerolag),
        
        if(isfield(S,'hmm_states')),
            sts=unique(S.hmm_states);
            for i=1:length(sts),
                tinds=(S.hmm_states==sts(i));
                dat1=S.dat(tri,tinds)';
                dat2=S.dat_seed(tri,tinds)';

                S.dat(tri,tinds)=dat1-dat2*(dat2\dat1);

                if(isfield(S,'dat_seed2')),
                    dat1=S.dat(tri,tinds)';
                    dat2=S.dat_seed2(tri,tinds)';

                    S.dat(tri,tinds)=dat1-dat2*(dat2\dat1);
                end
            end;     
        else,
            dat1=S.dat(tri,:)';
            dat2=S.dat_seed(tri,:)';

            S.dat(tri,:)=dat1-dat2*(dat2\dat1);
        end;        

    end;
    
    if(isfield(S,'tfdat_seed')),
        dat_seed=S.tfdat_seed;
    else,
        signal_h = hilbert(S.dat_seed(tri,:));

        % square the hilbert transform and take the
        % +ve square root to get the estimated envelope
        dat_seed=log(sqrt(signal_h.*conj(signal_h)))';
        dat_seed=dat_seed(time_indices);

        if(S.time_moving_av_win_size>0)
            fsample=1/(S.times_in(2)-S.times_in(1));        
            windowSize = S.time_moving_av_win_size*fsample;

            dat_seed_new=moving(dat_seed,round(windowSize));
            
            % figure;plot(dat_seed_new);ho;plot(dat_seed,'r');
            
            dat_seed=dat_seed_new;
            
        end;  

        
    end;
    
    signal_h = hilbert(S.dat(tri,:));

    % square the hilbert transform and take the
    % +ve square root to get the estimated envelope
    dat = log(sqrt(signal_h.*conj(signal_h)))';
  
    dat=dat(time_indices);
    if(S.time_moving_av_win_size>0)
        fsample=1/(S.times_in(2)-S.times_in(1));        
        windowSize = S.time_moving_av_win_size*fsample;

        dat=moving(dat,round(windowSize));
    end;
    

    cs=corrcoef(dat, dat_seed);cs=cs(1,2);
    
    conn2(tri,1)=cs;
    
    x=normalise(dat);
    y=normalise(dat_seed);

    pinvx=pinv(x);
    pe=pinvx*y;
    r=y-x*pe;
    vr=diag(r'*r/(size(y,1)-size(x,2)));
    vrp=pinv(x'*x)*vr;
    cs=pe/sqrt(vrp);
  
  
    
    conn(tri,1)=cs;        
           
end;

